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Application of the Probability Event Horizon filter to 
constrain the local rate density of binary black hole 
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ABSTRACT 

The temporal evolution of the gravitational wave background signal resulting from 
stellar-mass binary black hole (BBH) inspirals has a unique statistical signature. We 
describe the application of a new filter, based on the 'probability event horizon' (PEH) 
concept, that utilizes both the temporal and spatial source distribution to constrain 
the local rate density, ro, of BBH inspiral events in the nearby Universe. Assuming 
Advanced LIGO sensitivities and an upper rate of Galactic BBH inspirals of 30 Myr -1 , 
we simulate GW data and apply a fitting procedure to the PEH filtered data. To 
determine the accuracy of the PEH filter in constraining ro, a comparison is made 
with a fit to the brightness distribution of events. We apply both methods to a data 
stream containing a background of Gaussian distributed false alarms. We find that the 
brightness distribution yields lower standard errors, but is biased by the false alarms. 
In comparison the PEH method is less prone to errors resulting from false alarms but 
has a lower resolution as fewer events contribute to the data. Used in combination, the 
PEH and brightness distribution methods provide an improved estimate of the rate 
density. 
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1 INTRODUCTION 

The LIGO gravitational wave (GW) detectors are cur- 
rently taking data at design sensitivity, and embark- 
ing on long science runs. Promising GW sources po- 
tentially detectable by LIGO are coalescing binary sys- 
tems cont aining neutron stars (NSs) and/or blac k hole s 



(BHs)- s eefFlanaean fc Hughes! (| 19981 ) . iMiller et al l [|2004h . 



iThornd l|l995h . Because of the enormous GW luminosity 
~ 10" 3 c 5 /G ~ 10 23 L©, binary black hole (BBH) inspirals 
are among th e most promising candidates for a first detec- 
tion of GWs (iBaker et al.ll2002l: iFlanaean fc Hughes! 1 19981 : 



ISathvaprakashl 120041 : iBureer et al.ll200 



For LIGO-type detectors, even highly energetic BBH 
inspirals are predi cted to be detected at a rate of only 
(1CT 3 - 0.6)yr _1 (|Cutler and Thornd [2002] ') . However, the 
next generation of interferometric detectors, planned to 
go online in the next decade, should be sensitive to an 
abundance of sources, with event rates 1000 times greater 
than current detectors. These Advanced' interferometers 
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will provide a new window to the cosmos not accessible by 
conventional astronomy. 

The introduction of advanced detectors will allow us 
to detect BBH inspiral events out to a distance of z ~ 0.4. 
Events within this volume will contribute to the low 
probability 'popcorn' component of the astrophysical GW 
background for BBHs. It has been shown that the temporal 
and spatial distribution of transient sources which form this 
part of the GW background can be des cribed by the Prob- 
ability Event Horizon (PEH) concept of lCoward fc Burman! 
120051 ). 

The PEH describes the temporal evolution of the 
brightness of a class of transient events. The probability 
of a nearby event accumulates with observation time, so 
that the peak event amplitude has a statistical distribution 
dictated by the rate density and spacetime geometry. This 
feature provides a tool to model t he detectabilit y of a 
distribution of transient GW sources. ICoward et ail (120051 ) 
used the PEH to model the detectability of NS inspirals as a 
function of observation time assuming LIGO and Advanced 
LIGO sensitivities and reasonable estimates of the local 
rate density of events in the nearby Universe, ro. 
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In this study we consider the use of the PEH method 
to determine the local rate density of BBH coalescence 
events. We use a cosmological model to create synthetic 
data corresponding to four months data at advanced LIGO 
sensitivity. We then attempt to recover the assumed rate 
density. 

First we use the conventional approach of studying 
the brightness distribution of all events. This method does 
not make use of the time evolution. It simply utilizes the 
number-amplitude distribution and fits the observations to 
the rate density. We compare this with the PEH method, 
which utilizes the temporal distribution of events. For 
simplicity, we use a standard candle approximation to 
model our source population. We note that a network of 
detectors would allow us to exploit a special property of 
compact binary inspiral events, for which the chirp signal 
provides a measure of the luminosity distance to the source, 
enabling such sources to be treat e d as standard candles 
JSchutz)ll986l ; ISathvaprakashl 12004 IChernoff fc Finrj|2003t 
iFinnll 19931 ) 



To test both methods, we simulate a candidate popu- 
lation of BBH inspiral events by approximating the output 
data stream resulting from matched filtering. A simplified 
detection model (described in Section 5.1) is used which as- 
sumes Gaussian detector noise. We use the candidate pop- 
ulation to show that the temporal evolution of events has 
a unique statistical signature and exploit this signature to 
constrain ro- To determine the effect of detector efficiency 
we apply both methods to data corresponding to high and 
low false alarm rates. We also consider the effect of different 
star formation rate (SFR) models and determine bias intro- 
duced by different SFR evolution functions. 

The paper is organized as follows: In Section 2 we re- 
view event rate predictions and then discuss the evolution 
of the event rate of BBH inspirals in Section 3. In Section 4, 
we explain the PEH concept and show how it can be used to 
probe the source rate density and to differentiate between 
different astrophysical populations of sources. We describe 
the simulation of candidate events in Section 5 and use these 
data in Section 6 to determine an estimate of ro using the 
brightness distribution of sources. In Section 7 we present a 
method for extracting PEH data and use least-squares fit- 
ting to constrain ro. We present our results in Section 8 and 
in Section 9 summarize the key findings and discuss how this 
work can be extended. 



2 EVENT RATE ESTIMATIONS 

The local rate density of a particular astrophysical GW 
source is fundamental to estimating the number of poten- 
tially observable events. Usually defined within a volume 
spanning the Virgo cluster of galaxies, the local rate den- 
sity, ro, is determined using estimated source rates within a 
larger fixed volume of space. 

In the case of NS-NS inspirals, current rate estimates 
rely on a small sample of sources. The discovery of the 
double pulsar PSR J0737-3039, with estimated coalescence 
time ~ 87 Myr, increased the estimated inspiral rate for dou- 
ble NSs in our G alaxy by about an order-of-magnitude to 
20 - 300 Myr -1 l|Burgav et alj|2003l ; iKalogera etai1l2004 



with respect to earlier estimates l|Kalogera et all l200ll ; 
|Phinnevlll99ll ). 

The rates of BBH systems are even more uncertain. 
Because these systems have not been observed directly, 
their evolutionary parameters can be obtained only though 
population synthesis, which predicts Galactic event rates 
^ 10 -1 - 80 Myr _1 . However, GW emissions from BBH will 
be detectable out to much greater dis tances than other sys - 
tems of coalescing compact objects (|Sathva prakash 2004 . 
The rates of BH-NS inspirals are of a similar range to those 
of BH - BH systems, but have lower expected detection rates 
as a result of their less energetic emissions. 

In this study, we use the Galactic BBH coalescence 
rate TZ^ H ~ 30 Myr -1 , obtained from th e standard model 
in the population synthesis calculations of iBelczvnski et"aH 
(2002). We note that this rate is an upper limit. We convert 
this to a ra te per unit vo lume, ro, using the conversion factor 
10 -2 from lAndol (|l999h for the number density of galaxies 
in units of Mpc -3 , yielding the reference value of the local 
rate density, fo = 0.3 Myr -1 Mpc -3 , which we will employ 
in this study. 

Rate estimates will benefit greatly from the introduc- 
tion of advanced GW detectors, which will allow unobscured 
source counts to be conducted to almost cosmological vol- 
umes. A network of detectors will improve the sky coverage 
and source localization. A network may employ coherent 
analysis, in which synchronized detector outputs are merged 
l|Finnll200lh before a search for a common pattern, or alter- 
natively, a coincidence analysis, in which ind ividual events 
from different detectors are correlated in time l|Arnaud et al.l 
2002). However, as a result of the non-uniform antenna 
patterns, even a network of three detectors of similar sen- 
sitivity will have diffi culty obtaining maximum efficiency 
jArnaud et all l2003al ). In addition, the efficiency of a de- 
tector network for a particular source type will depend on 
the false alarm rate and the signal-to-noise threshold for de- 
tection - therefore high number counts may be balanced by 
an increased rate of false alarms. We will investigate this in 
Section 6. 



3 THE BBH COALESCENCE RATE 
EVOLUTION 

3.1 The event rate equation 

For standard Friedman cosmology a differential event rate 
in the redshift shell z to z + dz is given by: 



dVr e(z) 
dR — — : — — dz . 



(1) 



dz 1 + z 

where dV is the cosmology-dependent co-moving volume el- 
ement and R(z) is the all-sky (4-zr solid angle) event rate, 
as observed in our local frame, for sources out to redshift 
z. Source rate density evolution is accounted for by the di- 
mensionless evolution factor e(z), normalized to unity in our 
local intergalactic neighbourhood, and ro is the z = source 
rate density. The (1 + z) factor accounts for the time dila- 
tion of the observed rate by cosmic expansion, converting a 
source-count equation to an event rate equation. 

The cosmological volume element is obtained by ca lcu- 
lating the luminosity distance from (cf. iPeebles! (| 19931 ). p. 
332) 
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Figure 1. The solid line shows the dimensi onless SFR density 
evolut ion factor e(z) for the SFR model SF2 of Porci ani &; Madaul 
[|200fl ) in the flat-A (0.3, 0.7) cosmology. To allow for the average 
coalescence time of BBH systems, the dashed line shows the effect 
of a time delay of 1 Gyr on e(z). 



d L (Y) = (1+z)- 
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(2) 



(3) 



(4) 



#o J h(z') ' 
and using (Porciani and Madau 2001, eq. 3) 

aV _ 47re d h 2 (z) 
dz ~ H {l + z) 2 h(z) ' 

The normalized Hubble parameter, h(z), is given by 

h( z ) = h(z)/h = [n m (i + zf + n A ] 1/2 

for a 'fiat-A' cosmology (fi m + Qa = 1). We use Q m = 0.3 
and Qa = 0.7 for the z — density parameters, and take 
Ho = 70 km s _1 Mpc -1 for the Hubble parameter at the 
present epoch. 



3.2 The source rate evolution of BBH 

For the source rate evolution factor, e(z), we employ three 
star formation rate (SFR) models and a non-evolving SFR 
density for comparison. To simulate a candidate popula- 
tion of BBH inspiral ev ents we employ th e obse rvation- 
based SFR model SF2 of iPorciani fc Madaul (|200lh . Based 
on observed rest-frame ultraviolet and Hcv luminosity den- 
sities, this model includes an allowance for uncertainties in 
the amount of dust extinction at high z. In order to con- 
strain ro by least-squares fitting to the simulated data, we 
will use three additional models: a non-evolving SFR den- 
sity jno^el_o^tahied bysetting e(z) = 1, the model SF1 



sit y model obtained by setting e(z) = 1, the model at 1 
of iPorciani fc Madaul (200 1) and the model SH, based on 



an analytical fit to hydrodyna mic simulations conducted by 
ISpringel and Hernqu ist (2003) in a fiat-A cold dark matter 
cosmology; SF1 includes an upward correction for dust ex- 
tinction at high z. We re-scale SF1 and SF2, originally mod- 
elled in an Einstein-de-Sitter cosmology, to a fiat-A cosmol- 
ogy using the procedure outlined in the appendix of Porciani 
& Madau. 

ICoward et ail (|2005l ) assumed that the formation of 
double NS systems closely tracks the evolving star forma- 
tion rate. They based their assumptions on short merger 



redshift, z 

Figure 2. The all-sky BBH coalescence rate as a function of z 
using a Galactic rate TZq^ H ~ 30 Myr -1 and merger time of 1 
Gyr (Belczynski et al. 2002). We employ the observation-based 
star f ormation rate models SF1 and SF2 of IPorciani &; Madaul 
ll200ll) , a cons tant (non-evolving) model and a simulation-based 
model, SH, of lSpringel and Hernauid l|2003l l. 
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Figure 3. The all-sky rate of BBH coalescence as a function of z 
as previously shown in Figure[2] but within a distance of z = 0.4, 
corresponding to the Advanced LIGO detectability horizon for 
these sources. 



times and showed that a time delay of up to 5 Gyr had min- 
imal effect on the differential rate of events at low redshift. 
However, in comparison with NS-NS systems, for which the 
distribution of m erger times has a large range, peaking at 
around 0.3 Myr |Belczvnski et all 120021 '). the merger time 
distributions of BBH systems are predominantly skewed to- 
wards longer merger times - fr om about 100 Myr to the 
Hubble time (|Bulik et al.ll2004bl ). 

Belczynski et al. used population synthesis methods to 
calculate the properties and coalescence rates of compact 
binaries. They used a range of different scenarios defined by 
the initial physical parameters of the binary system such 
as component masses, orbital separations and eccentricities. 
They also included properties which affect the evolutionary 
channels of the system, such as mass transfer, mass losses 
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due to stellar winds and kick velocities. Their calculations 
implied that, compared to NS-NS systems, the wider orbits 
and stronger dependence of merger time on initial separation 
for BH-NS and BH-BH systems resulted in longer merger 
times, mostly in the range 0.1 to several Gyr (Belczynski et 
al. section 3.4.5). 

A standard model was defined by IBelczvnski et al] 

(2002) using a range of assumptions, including: a non- 
conservative mass transfer with half the mass lost by the 
donor lost by the system; a kick velocity distribution that 
accounts for the fact that many pulsars have velocities above 
500 km s -1 ; constant Galactic star formation for the last 
10 Gyr. For this model they found a distribution in BBH 
merger times that peaked at 1 Gyr. We therefore take this 
value as an average merger time and shift e(z) to reflect 
this delay tim e . Usin g the method described in section 2 of 
ICoward et all (|2005l ). we can convert e(z) to a function of 
cosmic time, t cos , using the relation 

Wz) = f [(l + i>(j')] 4 (i!'. (5) 
Jo 

We apply a 1 Gyr time shift and then convert back to a 
function of z. Figure Q] shows the factor e(z) for SF2 with 
and without the time delay. The result shows that within the 
range of advanced LIGO, z ~ 0.4, e(z) is not significantly al- 
tered. Therefore, although we include this time delay in our 
calculations, within the ranges of advanced LIGO detectors 
it will not have a significant influence on the results. 

Figure [5] shows the all-sky BBH coalescence rate, R(z), 
calculated by integrating the differential rate from the 
present epoch to redshift z. Using SF2, the rate continues 
increasing to distances well beyond the cosmological volume 
elements considered in this study. We therefore assume a 
universal rate of events, R v (usually defined as the asymp- 
totic value of the all-sky rate as z increases) of ~ 0.06 s - 
and a corresponding mean temporal interval r = 1/R V = 17 
s. 

Figure shows the all-sky BBH coalescence rate for 
sources within Advanced LIGO sensitivities. The potential 
detection horizon extends to z = 0.4 and the corresponding 
mean rate of events is ~ 2.5 x 10~ 4 s" 1 using SF2. If we 
assume all sources are composed of two IOMq black holes, 
this rate corresponds to around 874 events yr _1 with SNR 
> 8. 



4 THE PROBABILITY EVENT HORIZON FOR 
BBH COALESCENCE EVENTS 

4.1 The Probability Event Horizon 

The rate, as observed in our frame, of transient astrophysical 
events occurring throughout the Universe, is determined by 
their spatial distribution and evolutionary history. The dis- 
tribution of event observation times follows a Poisson distri- 
bution and the temporal separation between events follows 
an exponential distribution defined by the mean event rate. 
For cosmological events, the rate depends on the cosmology 
dependent volume and radial distance through redshift, z. 
We assume that an observer measures both a temporal loca- 
tion and a 'brightness' for each event, where the brightness 
is determined by the luminosity distance to the event. 

It follows from these assumptions that the probability 



for at least one event to occur in the volume bounded by z, 
during observation time T at a mean rate R(z) at constant 
probability e is given by the exponential distribution: 

p(n> l;R{z),T) = l-e- R( * )T = e, (6) 

1 — e~ RT being the probability of at least one event occurring 
(see Coward & Burman 2005). 

For equation ([6| to remain satisfied as observation time 
increases, the mean number of events in the sphere bounded 
by z, N e = R(z)T — |ln(l — e)|, must remain constant. The 
PEH is defined by the redshift bound, zf EH (T), required to 
satisfy this condition. 

We note that the for events occurring within a volume 
bounded by about a few Gpc, the PEH is well approximated 
using Euclidean geometry and takes on a simple analytical 
form with z now replaced by the radial distance in flat space: 

rf EH (T) = (3Af E /47rr ) 1/3 T- 1/3 , (7) 

where ro is a rate per unit volume. By setting e = 0.95, one 
can generate a threshold corresponding to a 95% probability 
of observing at least one event within zf EH (T), or alterna- 
tively, rf EH (T) for a Euclidean PEH model. We define a 
'null PEH', representing the 95% probability that no events 
will be observed within this threshold, by setting e = 0.05. 
When combined, we refer to these two PEHs as the 90% 
PEH band - that is 90% of events are expected to occur in 
the region enclosed by the two PEHs. By scaling some fidu- 
cial GW amplitude by d^(z), we express the 90% PEH band 
as 90% confidence bounds of peak GW amplitude against 
observation time. 

Figure|3]shows the Euclidean 95% PEH curves for BBH, 
using three different values of ro, differing by an order of 
magnitude. We note that for a particular transient GW 
source population, the PEH is intrinsically dependent on ro. 
The plot shows that for a particular source type, the PEH 
model can be used to estimate the value of ro- 



4.2 The PEH filter applied to astrophysical 
populations 

ICoward et al.l (|2005h outlined the PEH filter - a procedure 
used to extract a cosmological signature from a distribu- 
tion of events in redshift and time. The concept is very sim- 
ple: the longer one observes the greater the probability of a 
nearby event - the PEH filter quantifies this dependence as 
a means of probing the cosmological rate density. The PEH 
filter searches for the time dependence of event amplitudes 
which is imposed by their cosmological distribution. It is 
a non-linear filter applied to a body of data by recording 
successively closer events. We will apply the technique to a 
distribution of amplitudes, A, as a function of observation 
time, t, recording the successive (ti, Ai) that satisfy the con- 
dition Ai+i > Ai. The resulting events will be referred to as 
the 'PEH population'. 

By fitting to PEH data we can probe the local rate 
density of an astrophysical population. We note that just as 
brightness distributions can be used to separate and iden- 
tify different source populations, so also can different source 
populations be identified in PEH plots. 

Figure [5] demonstrates this property by comparing syn- 
thetic PEH populations of BBH and NS-NS inspirals. For 
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observation time (s) 

Figure 4. The GW amplitude 95% PEH curves for BBH using 
three values for the local rate density, ro. The PEH curve corre- 
sponding to the reference value of the local rate density, r o , used 
in this study is shown by the solid line. An increase in ro by an 
order of magnitude shifts the PEH curve towards earlier obser- 
vation times (dashed line). This implies an increased probability 
for detecting a local high-amplitude event. The dot-dashed line 
shows that a decrease in ro of the same magnitude has the reverse 
effect. 
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Figure 5. To illustrate that our simulated data is consistent with 
the PEH model, we show the 90% PEH thresholds (see Section 
4.1) and simulated PEH data for two different GW binary inspiral 
source types with varying event rates recorded within a redshift 
z. The PEH thresholds show the optimal SNR for an Advanced 
LIGO detector as a function of observation time. We assume a 
standard-candle approximation based on the fiducial distances for 
an optimal SNR of 8 for each source. The fiducial distances and 
universal event rates for the GW populations are: NS-NS inspirals 
at 200 Mpc and R u ~ 0.2 s _1 (represented by diamonds) and 
BBH at 2 Gpc and R u ~ 0.05 s _1 (indicated by triangles). 



the two source types we use a standard candle approxima- 
tion based on an optimal SNR of 8 assuming an Advanced 
LIGO detector with a fiducial distance of 200 Mpc for NS- 
NS inspirals and 2 Gpc for BBH inspirals. 

We show the 90% PEH bands for the two populations 
as a SNR. It is evident that if the PEH signature of an astro- 
physical GW background could be extracted from detector 
noise, the presence of different astrophysical backgrounds 
could be identified. Clearly this method will only work if 
the luminosities of the populations differ widely. In reality, 
different luminosities will also be associated with different 
waveforms, and there are likely to be more evident differ- 
ences than the simple luminosity effects. 



5 SIMULATING A CANDIDATE 

POPULATION OF BBH INSPIRALS 

5.1 The GW source model 

We will consider only the well understood inspiral stage of 
a BBH coalescence to provide quantitive data for our source 
model. To model the GW background from BBH inspirals 
it is sufficient to assume that the raw interferometer data 
is preprocessed by passing it through an optimal filter. This 
transforms e vents into approximate short duration Gaussian 
pulse signals (jAbbott et a l. 2004; S hawhan fc Ochsnerl2004h 
embedded in Gaussian detector noise. By injecting a popula- 
tion of simulated events into GW detector noise of Advanced 
LIGO sensitivity, we approximate the processed output data 
stream of a GW detector and then apply sub-optimal burst 
filtering to extract candidate events or 'triggers'. A more re- 
alistic detection pipeline will employ a range of templates 
representing different BBH mass configurations. To simplify 
our detection model, we assume a single BBH mass config- 
uration for our sources, representing the output of a fixed 
template. 

To represent the response to BBH events by optimal 
matched filtering, we adopt as a model waveform the sim- 
plified but robust fo rm used by lArnaud et all l|2003bl ) and 
lAbbott et~ai] ((2005) to model GW burst sources. This is a 
linearly polarized 5-ms duration Gaussian pulse, which ap- 
proximates to the form of the event triggers in processed 
data from the LIGO SI search for inspirals (jAbbott et al.l 
l2004l : IShawhan fc Ochsnerll2004l ). We note that we have cho- 
sen to ignore any additional secondary peaks which occur at 
high SNRs. 

The two polarizations of this signal will be given by: 



h+ (t) = A exp 



(t-to) 2 



2A Z 



h x (t) = 



(8) 



with amplitude A and half- width A; the time value of the 
signal maximum, to, is set to 10ms. The output response, 
h(t), will be a linear combination of the two polarizations 



h(t) = F+h+ +F x h > 



(9) 



where F + and F x are the antenna pattern functions, which 
are functions of sky direction, represented by the spherical 
polar angles 9 and 4>, and polar ization angle, tj), of the GW 
signals relative to the detector l|jaranowski et al]|l998l ). 

The filter response amplitude, A, will be dependent on 
the masses of the BBH system. For simplicity, rather than 
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• flat-A cosmology 

• source rate density 
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exponential distribution 
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Apply PEH 
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Figure 6. Flow chart outlining the simulation pipeline to gener- 
ate a candidate population of BBH inspiral events in interferom- 
eter detector noise. Individual inspiral events are scaled in am- 
plitude and time-dilated according to the random variable z, ob- 
tained from the probability distribution shown in Figure [7] The 
events are injected into simulated detector noise, with the tempo- 
ral separations between successive events following an exponen- 
tial distribution. Candidate signals are extracted using amplitude 
thresholding and robust sub-optimal filtering. 
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Figure 7. The probability distribution function for BBH inspirals 
throughout the Universe modelled using the same parameters as 
for the cumulative rate (see Fig.f3]l. For SF1, SF2 and the constant 
(non-evolving) model, the most probable events occur in z Ri 1 — 2. 
For SH, the most probable events occur significantly earlier, in 
z ?3 3 — 4, but with a flatter distribution. 



using a distribution of BBH masses, we use a standard source 
model. Using the BENCH^ code to model the detector noise 
spectrum, we approximate the GW amplitude for an opti- 
mum SNR of 80 at 200 Mpc. 

By assuming a standard BBH system we ignore the am- 
plitude distribution from the BH mass spectrum. However, 
this assumption is not a limitation to the PEH method be- 
cause the signature of an inspiralling binar y system con- 
tains a measu r e of it s luminosit y distance ([Schutz 19861 ; 
ISathvaprakash] |2004 iFinnl Il993l ; IChernoff fe Finn! I2003T I. 
Therefore, for a network of GW detectors, the analysis de- 
scribed in this paper could be repeated using luminosity 
distances rather than amplitudes. 



5.2 Simulation of GW interferometer data 

The simulation pipeline used to generate a cosmological GW 
population of BBH inspiral triggers in GW detector noise is 
shown in Figure [6] We use the BENCH code to calculate the 
noise sensitivi ty curve correspondin g to our source model. 

Following lArnaud et al.l(|2003bl ). we define the standard 
deviation of the detector noise as: 



a = /lrms\/ /o/2/ c 



(10) 



with /i rms the root-mean-square value of the advanced LIGO 
noise curve at rest frame frequency / c , and f the sampling 
frequency. We assume the noise is white Gaussian with zero 
mean. 

The amplitude and duration of each potential BBH in- 
spiral event is defined by the random variable z, generated 
from a probability density function P(z) for these events. 
We obtain P(z) by normalizing the differential event rate, 
dR/dz, by R u , the Universal rate of BBH inspiral events, 
integrated throughout the cosmos, as seen in our frame (cf. 
Coward & Burman 2005, section 3): 



P{z)dz = dR/R v 



(11) 



Equation 1111 defines P(z) dz as the probability that an 
observed event occurred in the redshift shell z to z + dz. In 
Figure[7]we present curves for the several star formation rate 
models. We see that the most probable events will occur at 
z ~ 1 — 2 for SF1 and SF2. The corresponding cumulative 
distribution function C(z), giving the probability of an event 
occurring in the redshift range to 2, is the normalized 
cumulative rate : 



C(z) = R(z)/R L 



(12) 



We use C(z) to simulate events and their associated 
redshifts and hence the GW amplitude of each injected can- 
didate event, h(t), is computed from equation Random 
values of the antenna response variables, 8, <f> and ip, are 
simulated and h(t) is scaled inversely by the luminosity dis- 
tance c(l(z). The signal duration is time-dilated by the factor 

(1 + 2). 

The temporal distribution of events in our frame is 
stochastic and the separation of events is described by Pois- 
son statistics. The time interval between successive events, 



1 The program BENCH can be obtained from the URL 
http://ilog.ligo-wa.caltech.edu: 7285/advligo/Bcnch 
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Figure 8. The top panel shows an optimally orientated burst 
signal (a 5- ms Gaussian pulse) at t fa 0.4 s with an optimal SNR 
of 8 when embedded in Gaussian noise corresponding to 200 Hz 
of the Advanced LIGO noise spectrum. The bottom panel shows 
the response of the mean filter in terms of a SNR. The mean filter 
response to this signal is about 70% of that of an optimal filter. 
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Figure 9. The effective SNR performance of three GW burst 
filters, as functions of z, in identifying the optimally orientated 
Gaussian pulse signal shown in Fig. \E\ Compared with the per- 
formance of an optimal filter in Gaussian detector noise with Ad- 
vanced LIGO sensitivity, the mean and norm filters provide useful 
performance. Noise samples are generated for each z - SNR fluc- 
tuations are the result of RMS variations between noise samples. 

r, will therefore follow an exponential distribution. Succes- 
sive waveforms, with amplitude and distance defined by the 
random variable z are generated and injected into a data 
pipeline at time intervals r. 

5.3 Candidate signal extraction 

We restrict ourselves to the well understood coalescence 
phase of BBH sources and as discussed earlier matched fil- 
tering provides processed data whereby candidate events can 
be approximated as Gaussian bursts embedded in a back- 
groun d of noise l|Abbott et al. 1 12004 IShawhan fc Ochsnerl 
2004). Candidate events will usually be selected on the 
basis of coincidence analysis and a \ 2 waveform consis- 



tency test performe d between template and filtered output 
jAbbott et al.ll2004l ). 

In our simulation pipeline we inject signals directly into 
Gaussian detector noise, thereby contaminating potential 
triggers; this makes a simple thresholding procedure insuffi- 
cient to extract good candidates. A combination of matched 
filtering and 'pulse' detecti on techniques such as 'burst fil- 
ters' has been suggested bv lPradier et al.l (|200ll ) as a means 
of increasing the final SNR for GW inspiral events. We there- 
fore employ burst filtering to extract a population of event 
candidates. 

For simplicity, we employ the mean filt e r, a hig hly ro- 
bust, linear filter developed bv lArnaud et al.l l|2003bl ). which 
operates in the time domain by calculating the mean of the 
data, Xi, in a sliding window of sample width N: 

1 j+JV-i 
i=j 

Figure [8] shows the response of this filter to an opti- 
mally orientated Gaussian pulse (see equation [8| at a dis- 
tance of z = 0.4, j ust within the expected detection limit for 
Adva nced LIGO ijSathvaprakashl |2004| ; ICutler and Thornel 
2002). The response is displayed as a maximum SNR - the 
maximum value of the ratio of mean filter output when a 
signal is p resent to the standard deviation in the absenc e 
of a signal l)Arnaud et al.ll2003bl : iFTanaean fc Hug hes 1998). 
The maximum response of about 5.5 is around 70% of that 
using optimal filtering - for which a SNR of about 8 was 
obtained. This mean filter response is typical of values ob- 
tained during testing. 

Figure [9] shows the comparative filter performances ob- 
tained for our pulse model at different values of z operating 
in white Gaussian noise, comparable in amplitude to that of 
the 200 Hz region of the Advanced LIGO noise curve. We 
compare the responses of the mean filter, norm filter and 
a simple band-pass filter, with that of a Wiener filter, with 
optimal SNR given by 

with Sh denoting the one-sided noise power spectral density. 
We see that the mean filter is the sub-optimal filter with the 
best average response. The fluctuations in the response are a 
result of generating different noise samples for each z. These 
results are in a greement with tests o n tim e domain filters 
carrie d out by lArnaud et alj l|2003bl ) and iBizouard et al.l 
(2003). The performance level, coupled with the robustness 
of this filter, make the mean filter an ideal choice for our 
candidate searches. 

The simulation pipeline for the generation of candidate 
events consists of the following steps: 

1. A 160-s data buffer, representing the output p(t) of 
a single optimal filter, is continuously populated with 
Gaussian detector noise at a sampling frequency of 1024 Hz 
— in an actual application, this represents a down-sampling 
of interferometer data from 2 14 Hz. 

2. Potential candidate events corresponding to different 
cosmological distances are injected into the data buffer at 
exponentially distributed intervals r as described in section 
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5.2. These intervals are recorded. 

3. The 160-s data segments are split into 80x2 s slices. 

4. If a slice contains no fluctuations above a threshold of 
a — 3.2, it is rejected. Otherwise, a mean filter is applied to 
the slice - an event with maximum SNR > 4.2 is recorded 
as a candidate BBH inspiral event, and is added to a 
candidate event population, E. 

5. The PEH algorithm is applied to the candidate events 
to extract the consecutive running maximum amplitudes - 
these events we describe as the candidate PEH population, 
C. 

Both event triggers and injected signals are recorded for 
later analysis. 



5.4 The candidate event population 

Figure QJj] shows the injected events representing the GW 
background of BBH within z ~ 5 for four months of observa- 
tion time, corresponding to around 500,000 events — about 
2200 of these events were within z — 0.4, the detectabilty 
horizon of Advanced LIGO, a value which is within the up- 
per limit of 8000 yr" 1 set bv lBelczvnski et all l|2002l) . 

Events from the peak of the probability distribution 
function (see Fig. dominate short observation times. As 
observation time increases, the rarer events at both high 
and low redshift become more numerous. The high-z events 
will be buried in detector noise, and only potentially de- 
tectable by cros s-correlation between co-located detectors 
jMaggiord l200Ch . The background component from more 
luminous events, z < 0.4, could be detected by Advanced 
LIGO as individual burst events. When we apply the PEH 
algorithm to these types of data, the PEH population, C, 
will consist mostly of the rarer events at the top edge of Fig. 

ED 

Figure [TT] displays the output of our simulation pipeline 
for an observation time of four months, constituting 37,387 
events. Of this population, 2173 candidates were identified 
as injected signals - the remainder are false alarms. At early 
observation times noise events dominate and the false alarm 
rate is high. As observation time increases, a greater pro- 
portion of the candidate signals are GW burst events. 



■.- 
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observation time (days) 

Figure 10. The simulated BBH inspiral population assuming 
a local rate density of 0.3 Myr - 1 Mpc - 3 for 4 months of data, 
corresponding to 527,680 events. Note how the GW amplitudes 
incorporate both higher and lower magnitude extremes as obser- 
vation time increases - corresponding to smaller and larger z. 
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Figure 11. The candidate population of BBH for 4 months of 
observation time - around 37,000 events — represented by grey 
squares. Of these candidates, the 2173 events that correspond to 
injected signals are shown as black diamonds. The remainder are 
classed as false alarms. 



5.5 The PEH population 

Figure [12] shows the PEH population of BBH events ex- 
tracted from one 4- month data segment. Each point repre- 
sents the maximum amplitude of events observed as a func- 
tion of observation time. Also shown is the BBH inspiral 
amplitude PEH of Fig. [S] and the Gaussian noise amplitude 
PEH threshold, which models the progressive maximum am- 
plitude growth of the Gaussian detector noise used in our 
simulation (see section 4.2). The small gradient of the noise 
PEH in comparison with the astrophysical PEH curve high- 
lights the fact that the temporal evolution of the Gaussian 
detector noise, highly dependent on the sampling frequency, 
is much slower, a result o f the low probability of events in the 
tail of the distribution (jCoward et al.ll2005T ). Similarly, the 



narrower 90% threshold for the Gaussian-noise PEH is also 
a result of tighter constraints on the amplitude distribution 
for the Gaussian noise model. 

For early observation times (< 30,000 s) the PEH popu- 
lation is dominated by false alarms, lying close to the Gaus- 
sian noise PEH. As observation time increases, the astro- 
physical PEH population begins to dominate. An additional 
threshold, shown by the dot-dashed line, is constructed by 
combining the null-PEH curves from both populations. This 
threshold is the 95% amplitude upper limit for our BBH 
inspiral PEH population - this curve represents the result 
of obtaining maximum amplitudes for both the noise and 
sources. 
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Figure 12. The PEH population of BBH for 4 months of obser- 
vation time at advanced LIGO sensitivity are shown by crosses, 
with the solid lines representing the associated 90% amplitude 
PEH thresholds (see Fig[5]l. In addition, we show the 90% Gaus- 
sian noise amplitude PEH curves (dashed lines), which illustrate 
the temporal amplitude evolution of Gaussian detector noise. The 
dot-dashed line represents the absolute amplitude upper limit at 
95% confidence for signal + noise. 

6 THE LOG N - LOG A DISTRIBUTION 

We now want to determine the key astrophysical parame- 
ter, the source local rate density ro, from data such as that 
presented in the previous section. First we will consider the 
conventional method of determining source rate densities 
based on the brightness distribution of sources. This will 
allow us to quantify the effectiveness of the PEH in deter- 
mining the local rate density o f BBH events . The log N- 
log P source count d istribution l|Guetta et al. I l2005l ; iTotanl 
1 19991 ; ISchmiddl200lh is the number of events N(> P), of 
luminosity L and peak flux P, within a maximum redshift, 
Zmax(-£/, -Piim) 5 recorded by a detector of flux limit Pu m . Such 
a distribution can be fitted to a predicted log N-log P curve 
to obtain rate estimates or to constrain the luminosity func- 
tion. 

For our standard-candle GW population, we convert 
peak flux to a maximum GW amplitude, A, yielding a log 
N-log A distribution of the form: 

rz(L,A) dR 

N(>A)= / ro^dz (15) 

Jo 

where to represents the observation time and the differen- 
tial event rate, d_R/dz, is given by equation[T] In comparison 
with the log N-log P distribution, which has a gradient of 
—3/2 under a Euclidean geometry, the log N-log A has a 
slope of —3. The curve includes a noise component which 
approximates the average contribution from detector noise 
and is scaled by a factor 0.4, the mean valu e of the ant enna 
response function for a single GW detector l|Finnl ll993l). By 
fitting to 100 synthetic data sets, we find these approxi- 
mations introduce a systematic error of ±6% to the final 
estimates. 

To fit the candidate event population, E, against the log 
N-log A curve, we consider two scenarios: firstly an idealized 
case, in which the detector has correctly resolved all the 



Data stream 


Q T?T3 , , , ,,-1^1 

orK model 


Estimate of ro 






(in units of fo) 


i 


SF2 


1.00 ±0.07 


i 


constant 


1.74 ±0.11 


i 


SF1 


0.99 ± 0.07 


i 


SH 


1.46 ± 0.10 


2 


SF2 


1.01 ± 0.07 


2 


constant 


1.67 ±0.10 


2 


SF1 


0.99 ± 0.07 


2 


SH 


1.48 ±0.10 



Table 1. The results of least squares fitting to the log N - log A 
distributions of two independent data steams, each representing 4 
months of observation time. The data consists of 2273 events for 
data steam 1 and 2173 events for data stream 2. The results for 
data stream 1 are shown in Fig |13l The data streams represent the 
output of a perfect detector - all false alarms have been dismissed. 
The SFR model is the one used in the fit; the estimated ranges 
of ro are given at 90% confidence. 



Data stream SFR model Estimate of ro 

(in units of fo) 



1 


SF2 


1.36 ± 0.08 


1 


constant 


1.96 ± 0.13 


1 


SF1 


1.35 ± 0.08 


1 


SH 


1.80 ± 0.12 


2 


SF2 


1.28 ± 0.08 


2 


constant 


1.86 ± 0.13 


2 


SF1 


1.27 ±0.08 


2 


SH 


1.77 ±0.13 



Table 2. The same as for Table 1, but for data with a high false 
alarm rate. The data consists of 38241 events for data stream 1 
and 37873 events for data stream 2. The results for data stream 
1 are shown in Fig 1141 

injected events in E; secondly a suboptimum case, in which 
the data consists of both injected signals and false alarms. 
To implement the first scenario we eliminate any false alarms 
by including only event triggers that correspond to injected 
signals. For the second scenario we use the whole candidate 
population of event triggers and false alarms. In section 8 
we will apply the PEH method to both scenarios, thereby 
allowing a direct comparison. 

Figure [13] shows the results of fitting against an ideal- 
ized data set. For data streams 1 and 2 we obtain estimates 
of (1.00±0.07)f and (1.01 ± 0.07)f within 90% confidence 
using a non-linear regression model based on SF2, the model 
used to simulate our data. Table 1 shows that the estimates 
obtained using fits based on four different SFR models re- 
cover the true rate to within an average of 46%. These es- 
timates represent the upper limits obtainable by applying 
the log N - log A method to our simulated data set. We 
note that the choice of SFR model introduces a bias in our 
estimates of ro. These biases result from the SFR model de- 
pendence of the all-sky event rate, R{z), as illustrated in 
Figure [3] 

Figure [14] highlights the effect of false alarms on the 
data sample. At low amplitude, the distribution is aug- 
mented by the noise transients resulting in overestimates 
of the local rate density. We therefore use a thresholding 
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Figure 13. A least squares fit to simulated linear-linear GW 
data, using predicted brightness distributions based on different 
SFR models, with a single free parameter, rg. The data are candi- 
date events from an event population that correspond to injected 
signals — this represents an idealized data set, in which we as- 
sume all false alarms have been dismissed by the detector. The 
top panel shows the result of fits based on SF2 and a constant 
(non-evolving) model; the bottom panel shows fits based on SF1 
and SH. The estimates obtained using this fitting procedure are 
given in Table 1. 

procedure to dismiss all events below log A = -21.55. Us- 
ing a fit based on the SF2 model, we obtain estimates of 
(1.36 ± 0.08)f using data stream 1 and (1.28 ± 0.08)f for 
data stream 2. The estimates obtained using fits based on 
the four different SFR models are shown in Table 2. The 
consistency of each pair shows that the statistical errors are 
small. In addition to the biases resulting from the differ- 
ent SFR models used, for this suboptimum case there is a 
significant bias due to the effect of false alarms in the distri- 
bution. In comparison, one might expect the PEH method, 
which considers only the most energetic events as a function 
of time, should be less sensitive to false alarms due to noise, 
but may have less precision since not all events are utilized. 
This is considered in the analysis below. 



7 FITTING TO THE PEH POPULATION 

The dependence of the PEH on ro, as demonstrated in Fig. 
21 provides a means of fitting the PEH curves to a candidate 
PEH population, C, and estimating ro. However, an obstacle 
to this procedure is highlighted in Fig. 1121 which shows that 
C is dominated by detector noise at early observation times 
(t < 30,000 s). This reduces the samples available to fit to 
the amplitude PEH curves. 

To reduce the inclusion of false alarms in the fitting pro- 
cedure, we only include the sample of C above 30,000 s. This 
threshold corresponds to the change in gradient of the PEH 
distribution as injected events become dominant over the 
Gaussian noise components at long observation times (see 
Fig. 11 2 pi . This choice of threshold will vary for different as- 
trophysical populations, depending on the emission energies 
and rates. After applying this threshold, we expect a total 
PEH population, C, of about 7-11 events for four months of 
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Figure 14. The non-linear regression models of Fig. 1131 are fitted 
to the full candidate population, consisting of both injected events 
and false alarms. We fit to the data above a threshold of log A 
= -21.55. The estimates obtained using this fitting procedure are 
given in Table 2. 



Test Subset Sample space of KS Probability 



combination test samples P(> Zn) 



1 1x4 month 8, 8 0.36 

2 2x2 month 10, 11 0.59 

3 3 x 1.3 month 17, 18 0.74 

4 4x1 month 22, 21 0.81 



Table 3. The results obtained from 2D Kolmogorov-Smirnov 
tests between two 4-month samples of BBH data, extracted from 
the same simulated data stream. Each sample is split into dif- 
ferent numbers of subsets of equal duration and recombined to 
increase the sample size. We see that the KS statistic improves 
with sample space for these combinations. 



observation time. We can compensate for any loss of events 
by this thresholding procedure, by utilizing the signature of 
the PEH to increase our sample space as discussed below. 

We can increase the sample space of C by splitting our 
candidate event population, E, into I subsets of equal obser- 
vation time. By applying the PEH filter to each subset and 
recombining, we form a more highly populated sample C - 
the temporal duration of which will now correspond to the 
length of each individual subset. 

The available subset configurations are given by Si, i — 
l,...,l. As the PEH is independent of when the detector 
is switched on, the PEH signature imprinted within each 
subset, Si, will be set by the subset's duration. When we 
combine all Si, to form C, the individual signatures will 
be subsequently imprinted within the overall sample space. 
This useful procedure enables us to increase the statistical 
sample. 

When splitting and recombining C, there is a deli- 
cate balance between improving the statistics by increas- 
ing the sample space, and reducing it by decreasing the 
duration. The most efficient way to increase the popula- 
tion of C, whilst retaining the embedded statistical signa- 
ture of the PEH can be deter mined using a two-dimensiona l 
Kolmogorov-Smirnov test (jFasano fc Franceschinil Il987l ). 
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Figure 15. A least squares fit to simulated GW data using a 36% 
PEH curve based on SF2 as a non-linear regression model. The 
data correspond to an idealized situation in which all false alarms 
have been dismissed. Using a PEH fit based on SF2 we obtain an 
estimate of (1.26 ± 0.56)fo for data stream 1 and (0.84 ± 0.53)ro 
for the data stream 2. As an additional test of the PEH model, 
the dashed lines show that the 90% PEH thresholds are a good 
fit to the data. 



This test determines the probability of two sets being drawn 
from the same population distribution. 

We calculated the two-dimensional Kolmogorov- 
Smirnov probabilities P KS obtained using two consecutive 4- 
month samples taken from the same simulated data stream. 
A high value of P KS will provide evidence of the statisti- 
cal compatibility of two samples. Table 3 shows that a KS 
probability of 81% was obtained by splitting the data into 
four 1-month subsets and recombining. We therefore use this 
configuration to produce C, the PEH data to which we will 
apply our fitting procedure. 

It should be noted however, that the best choice 
of I, corresponding to the sample configuration with the 
strongest statistical PEH signature, may vary for different 
candidate event populations, E, depending on the observa- 
tion time. For example, in performing the same test on three 
months of data we found that a larger P KS was obtained us- 
ing I = 3 than for I = 4 - the result of a loss in the PEH 
signature for samples of shorter temporal duration. For this 
paper, to illustrate the PEH technique and avoid any ad- 
ditional complications, we maximize the sample space by 
using a 4-month sample of data. 

To constrain the local rate density we apply a least- 
squares fit to the candidate population C using an amplitude 
PEH curve as a linear regression model with free parameter 
ro. When equation [6] is fitted to the data it is necessary 
to determine the value of e that correctly estimates ro- This 
value, equal to 0.36, was experimentally verified by fitting to 
1000 synthetic data sets using e as the only free parameter. 
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Figure 16. The same fitting function as used in Figure IT51 is ap- 
plied to data which include both injected events and false alarms. 
We obtain an estimate of (1.43 ± 0.61)frj for data stream 1 and 
(0.83 ± 0.51)r for the data stream 2. We again fit the 90% PEH 
thresholds to the data, shown by the dashed lines. 



Data stream 


SFR model 


Events 


Estimate of ro 








(in units of ?o) 


1 


SF2 


19 


1.26 ±0.56 


1 


constant 


19 


1.43 ±0.63 


1 


SF1 


19 


1.25 ±0.56 


1 


SH 


19 


1.37 ±0.60 


2 


SF2 


18 


0.84 ± 0.53 


2 


constant 


18 


0.94 ±0.57 


2 


SF1 


18 


0.83 ±0.52 


2 


SH 


18 


0.92 ±0.56 



Table 4. The results of least squares fitting to the PEH distribu- 
tions of two independent data steams, each representing 4 months 
of observation time. The data represents that of an idealized situ- 
ation in which all false alarms have been dismissed. A PEH fit to 
data stream 1 is shown in Fig |15l The SFR model is the one used 
in the fitting function. The numbers of events used in each fit are 
shown along with the estimates of ro given at 90% confidence. 



Data stream 


SFR model 


Events 


Estimate of ro 








(in units of ro) 


1 


SF2 


21 


1.43 ±0.61 


1 


constant 


21 


1.61 ±0.71 


1 


SF1 


21 


1.42 ±0.63 


1 


SH 


21 


1.55 ±0.69 


2 


SF2 


19 


0.83 ±0.51 


2 


constant 


19 


0.94 ±0.55 


2 


SF1 


19 


0.82 ±0.51 


2 


SH 


19 


0.91 ±0.54 



Table 5. The same as for Table 1, but for data with a high false 
alarm rate. A PEH fit to data stream 1 is shown in Fig 1161 
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8 RESULTS 

In this section we present the results obtained by least- 
squares fitting to candidate PEH populations C using the 
36% amplitude PEH curve as a linear regression model with 
free parameter ro. We use the same two sets of data as in 
section 5 and consider again two different scenarios: firstly, 
an idealized case in which the detector operates at high ef- 
ficiency and has correctly dismissed all false alarms, and 
secondly, a suboptimum case in which data include false 
alarms. 

Figure [15] shows the results of a non-linear fit to both 
data streams for the idealized case in which all false alarms 
have been dismissed. This fit uses a PEH curve based on 
SF2, the model used to generate the data, as a non-linear 
regression model and yields estimates of (1.26 ± 0.56)fo for 
data stream 1 and (0.84 ± 0.53)fo for data stream 2. As 
an additional test of the PEH model, this figure shows that 
the data is well constrained by the 90% amplitude PEH 
thresholds. 

Estimates of ro obtained using PEH curves based on all 
SFR models are shown in Table 4. We see that we obtain 
mean estimates to within around a factor of 1.5 for data 
streams 1 and 2. As a result of the smaller data set used 
in the PEH fitting procedure, estimates can be sensitive to 
the distribution of data. This effect is highlighted in Figure 
[15] which shows that the PEH distribution of data stream 1 
tends towards the upper PEH threshold, producing higher 
estimates of ro. 

The uncertainties in these estimates are greater than 
those obtained using the brightness distribution. For exam- 
ple, using a fit based on SF2, the brightness distribution 
recovered the true rate within 7% and 8% for data streams 
1 and 2, while the equivalent estimates obtained using the 
PEH fit were within 82% and 69%. We note however, that 
the results obtained using the PEH method are not as promi- 
nently effected by bias due to different SFR models as those 
obtained using a brightness distribution. 

In Figure \W\ we see the result of a least-squares fit 
to data which contains both the injected events and false 
alarms. The numerical estimates determined by applying 
the PEH fitting procedure to these data are presented in 
Table 5 for both data streams 1 and 2. The mean estimates 
for both data streams are within a factor of 1.5 of the true 
value of ro- 

To determine the effect of false alarms on our estimates 
it is useful to eliminate the effects of using different SFR 
models in our fitting procedures. This can be best achieved 
by looking at the results obtained using SF2, the model used 
to produce the candidate data streams. We see that a PEH 
fit based on SF2 yields estimates of (1.43 ± 0.61)fo for data 
stream 1 and (0.83 ± 0.51)fo for data stream 2. We see that 
the inclusion of false alarms has resulted in a 24% and a 
3% decrease in accuracy for data streams 1 and 2 respec- 
tively. In comparison, estimates using an SF2 model to fit 
to the brightness distribution degraded by up to 44% for 
data stream 1 and 35% for data stream 2. 

These results imply that the PEH method is less prone 
to errors due to the inclusion of false alarms. This outcome 
arises because the PEH distribution is composed of the most 
energetic events as a function of observation time, thereby 
omitting most false alarms. We note however that the PEH 
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Estimate of ro 
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SF2 


1.00 ±0.06 
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constant 


1.73 ± 0.10 
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SF1 


0.99 ± 0.06 
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SH 


1.45 ± 0.09 


2 


SF2 


1.00 ±0.06 


2 


constant 


1.64 ±0.09 


2 


SF1 


0.98 ± 0.06 


2 


SH 


1.46 ±0.09 



Table 6. The constraints on the true local rate density, ro, ob- 
tained by combining the results of least squares fitting to the 
log N — log A distributions (Table 1) and the PEH distributions 
(Table 4) of the two independent data steams. The data sets rep- 
resent the output of a perfect detector — all false alarms have been 
dismissed. As previously, the SFR model is the one used in the 
fit; the estimated ranges of ro are given at 90% confidence. 



Data stream 


SFR model 


Estimate of ro 






(in units of fo) 


1 


SF2 


1.36 ± 0.07 


1 


constant 


1.94 ±0.12 


1 


SF1 


1.35 ± 0.07 


1 


SH 


1.79 ± 0.11 


2 


SF2 


1.26 ± 0.07 


2 


constant 


1.81 ± 0.12 


2 


SF1 


1.25 ± 0.07 


2 


SH 


1.72 ± 0.12 



Table 7. The same as for Table 6, but this time we combine the 
estimates of Tables 2 and 5 for data with a high false alarm rate. 

method has a lower resolution than that of the brightness 
distribution - a direct result of the smaller data sets used 
in the fitting procedure. This obvious disadvantage will be 
discussed in the next section. 

Tables 6 and 7 show the estimates obtained by combin- 
ing brightness distribution and PEH data. When compared 
to the results obtained using the brightness distribution in 
Tables 1 and 2, we see that including the PEH data improves 
the estimates by at least 1-2% in most cases. In addition, 
these results indicate that the PEH method can be employed 
as an additional test of consistency. 



9 CONCLUSIONS 

We have shown that the PEH filter allows an independent 
estimate of the rate density of BBH coalescence events de- 
tected by advanced GW detectors. The main results and 
their limitations are summarized below: 

(i) For a candidate population of BBH with Galactic in- 
spiral rate !ZfJi H ~ 30 Myr -1 , a fit to 4 months of inter- 
ferometer data was sufficient to obtain estimates of ro to 
within a factor of 2 at the 90% confidence level. 

By applying both brightness distribution and PEH 
methods to data streams with a high false alarm rate we 
find that the brightness distribution method is the more 
accurate if the detector is operating at high efficiency. For 
the case in which the data contains a large proportion of 
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false alarms, we find the PEH fit gives less bias. This is 
because this method naturally suppresses low-amplitude 
false-alarms. Of the two methods, the PEH method has 
lower resolution due to the fact that fewer events con- 
tribute to the data. The overall performance of the PEH 
filter suggests that it is accurate enough to be considered 
as an additional tool to determine event rate densities, 
particularly if the detector is operating at low efficiency. 
A combination of the PEH and brightness distribution 
methods provides two independent estimates of ro- The 
combination of both estimates provides a self consistent test 
and also increases the overall precision of the rate estimates. 

(ii) One disadvantage of the PEH filter, in comparison 
with fitting to the amplitude distribution, is that it uses 
only a small sample of the overall data set. We are investi- 
gating techniques to increase the sample size. Initial results 
suggest that applying the PEH filter in both temporal 
directions can increase the overall PEH population and 
improve the resolution and accuracy of the estimates. This 
technique will be particularly useful in data sets in which a 
large event occurs after a comparatively short observation 
time. 

(iii) We note that the value for reference local rate 
density, f o, obtained from t he po pulation synthesis calcu- 
lations of iBelczvnski et al l (|2002l ). is at the upper end of 
predictions. We plan to investigate the performance of the 
PEH method for lower values of fo for which we expect a 
smaller number of events in the PEH distribution. Initial 
calculations suggest that an order of magnitude decrease 
in fo will result in a PEH population of around 7-11 
events for a 4-month data set rather than the 15 - 22 events 
expected for f q used in this study. As discussed previously, 
techniques in which we can increase the PEH sample will 
be of great importance for astrophysical populations with 
lower rate densities. 

(iv) For the noise levels considered in this study, the 
PEH fits are only weakly affected by the inclusion of false 
alarms, but estimates using a brightness distribution are 
shown to degrade. This implies that the PEH method may 
be most effective when applied to data output from signals 
that are not well modelled, such as transient burst sources, 
for which we expect a substantial number of false alarms 
to be present. Recent developments in the modelling of 
GW emissions from core-collapse supernova (SNe) suggest 
that the detectable background population of such sources 
may be numero us enough to apply the PEH technique 
l|Ott et al.ll2006h . 

(v) The PEH filter, by definition includes the temporal 
distribution of the sources, hence provides a means of 
predicting the energies of future events. 

(vi) Our results show that at the detector sensitivity as- 
sumed here, the value of ro cannot be separated from evolu- 
tion effects, so that different SFR curves create bias in the 
estimated value of ro . For larger spans of data and m ore sen- 
sitive detectors such as EURO <|Sathvaprakashll2004h it may 
be possible to fit the PEH data or the brightness distribution 
to determine the entire function R(z). 
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